classdef Assembler < handle
    % AssemblyManager 有限元系统组装管理器
    % 负责协调Entity、Material和Element之间的数据流，组装全局系统刚度,质量矩阵,荷载向量
    
    properties
        % --- 系统组件 ---
        % 需要注入的信息
        ID
        Entities            % Entity对象数组
        Materials           % Material对象数组
        GlobalNodes         % 全局节点坐标表 [nNodes × dim]
        LoadSet             % 荷载对象
        
        % --- 全局系统矩阵 ---
        % 输出信息
        GlobalStiffness     % 全局刚度矩阵
        GlobalMass          % 全局质量矩阵 
        GlobalLoad          % 全局荷载向量
        
        % --- 系统信息 ---
        NumTotalNodes       % 总节点数
        NumTotalDofs        % 总自由度数
        
        % --- 日志信息 ---
        ElementCount
    end
    
    properties (Dependent)
        SystemInfo          % 系统基本信息
    end
    
    methods
        %% ===================== 构造函数 =====================
        function obj = Assembler(globalNodes, ID)
            % 构造函数
            % globalNodes: 全局节点坐标表
            % ID: 全局ID矩阵
            
            obj.GlobalNodes = globalNodes;
            obj.ID = ID;

            obj.NumTotalNodes = size(globalNodes, 1);
            obj.NumTotalDofs = max(ID(:));
            
            % 初始化系统矩阵
            obj.GlobalStiffness = sparse(obj.NumTotalDofs, obj.NumTotalDofs);
            obj.GlobalMass = sparse(obj.NumTotalDofs, obj.NumTotalDofs);
            obj.GlobalLoad = zeros(obj.NumTotalDofs, 1);
            
            obj.Entities = [];
            obj.Materials = [];
        end
        
        %% ===================== 系统组件注册 =====================
        function registerEntities(obj, entities)
            % 注册Entity对象数组
            validateattributes(entities, {'Entity'}, {'vector'});
            obj.Entities = entities;
        end
        
        function registerMaterials(obj, materials)
            % 注册Material对象数组
            validateattributes(materials, {'Material'}, {'vector'});
            obj.Materials = materials;
        end

        function registerLoadSet(obj, LoadSet)
            % 注册Entity对象数组
            validateattributes(LoadSet, {'LoadSet'}, {'vector'});
            obj.LoadSet = LoadSet;
        end
        
        %% ===================== 核心组装方法 =====================
        function assembleSystem(obj)
            % 组装完整的有限元系统
            
            fprintf('开始组装有限元系统...\n');
            
            % 1. 清空系统矩阵
            obj.clearSystemMatrices();
            
            % 2. 遍历所有Entity，创建并组装单元
            totalElements = 0;
            for i = 1:length(obj.Entities)
                entity = obj.Entities(i);
                entityElements = obj.assembleEntityElements(entity);
                totalElements = totalElements + entityElements;
            end
            
            fprintf('系统组装完成: %d 个单元, %d 个节点, %d 个自由度\n', ...
                totalElements, obj.NumTotalNodes, obj.NumTotalDofs);

            % 3. 组装荷载向量
            obj.assembleNodeLoadsToF()
        end
        
        function elementCount = assembleEntityElements(obj, entity)
            % 组装单个Entity中的所有单元
            obj.ElementCount = 0;
            
            entity.displayInfo();
            
            % 使用Entity的forEachElement方法遍历所有单元
            entity.forEachElement(obj.GlobalNodes, ...
                @(elementType, connectivity, materialIndex, elementCoords) ...
                obj.processSingleElement(elementType, connectivity, materialIndex, elementCoords));
            
            elementCount = obj.ElementCount;
            fprintf('  -> 实体 %d 组装完成: %d 个单元\n', entity.ID, elementCount);
        end
        
        function processSingleElement(obj, elementType, connectivity, materialIndex, elementCoords)
            % 处理单个单元的组装过程
            
            obj.ElementCount = obj.ElementCount + 1;
            
            try
                % 1. 获取材料模型
                material = obj.Materials(materialIndex);
                if isempty(material)
                    warning('未找到材料ID: %d，跳过单元 %d', materialIndex, elementCount);
                    return;
                end
  
                % 2. 创建对应的Element对象，调用单元工厂
                element = ElementFactory.createElement(elementType, connectivity, material);
                
                if isempty(element)
                    warning('无法创建单元类型，跳过单元 %d', obj.elementCount);
                    return;
                end

                % 3. 设置节点坐标
                element.setNodeCoordinates(elementCoords);

                % 4. 计算单元矩阵并组装到全局系统
                obj.assembleElementMatrices(element);
                
            catch ME
                warning('输出：单元 %d 组装失败: %s', obj.ElementCount, ME.message);
            end
        end
        
        function assembleElementMatrices(obj, element)
            % 组装单个单元的矩阵到全局系统

            % 获取单元自由度索引
            dofIndices = element.calculateElementLM(obj.ID);

            % --- 过滤掉被约束的自由度 ---
            % 找出哪些自由度是活动的（ID不为0）
            activeDofs = dofIndices > 0;

            % 只保留活动自由度的索引
            activeDofIndices = dofIndices(activeDofs);

            % 如果该单元所有自由度都被约束，则跳过组装
            if ~isempty(activeDofIndices)

            % 计算单元刚度矩阵 Ke
            % Ke 的维度是 [nDofsPerElement x nDofsPerElement]
            Ke = element.computeStiffnessMatrix();

            % --- 从Ke中移除被约束自由度对应的行和列 ---
            % 这保证了 Ke_filtered 和 activeDofIndices 的维度完全匹配
            Ke_filtered = Ke(activeDofs, activeDofs);

            % 组装到全局刚度矩阵
            % 现在索引都是有效的正整数，可以安全地进行组装
            obj.GlobalStiffness(activeDofIndices, activeDofIndices) = ...
                obj.GlobalStiffness(activeDofIndices, activeDofIndices) + Ke_filtered;
            end
%             % 4. 可选：组装质量矩阵
%             try
%                 Me = element.computeMassMatrix();
%                 obj.GlobalMass(dofIndices, dofIndices) = ...
%                     obj.GlobalMass(dofIndices, dofIndices) + Me;
%             catch
%                 % 某些单元可能不支持质量矩阵计算
%             end
        end

         %% ===================== 荷载组装 =====================
        function assembleNodeLoadsToF(obj)
            % LoadData: LoadSet 对象
            % ID: global ID matrix [nNodes x nDofsPerNode] (0 means constrained)

            for k = 1:length(obj.LoadSet.NodeLoads)
                entry = obj.LoadSet.NodeLoads(k);
                nodeID = entry.nodeID;
                dof = entry.dof;
                value = entry.value;

                if nodeID <= 0 || nodeID > size(obj.ID,1)
                    warning('assembleNodeLoadsToF: invalid nodeID %d', nodeID);
                    continue;
                end
                if dof <= 0 || dof > size(obj.ID,2)
                    warning('assembleNodeLoadsToF: invalid dof %d for node %d', dof, nodeID);
                    continue;
                end

                gDof = obj.ID(nodeID, dof);
                if gDof ~= 0
                    obj.GlobalLoad(gDof) = obj.GlobalLoad(gDof) + value;
                else
                    % 自由度被约束（ID=0），外力作用在约束自由度上通常
                    % 要么视为已知反力的一部分，要么在求解前警告
                    % 这里我们只记录警告，方便用户检查
                    warning('assembleNodeLoadsToF: applied load on constrained DOF node %d dof %d (ignored in F).', nodeID, dof);
                end
            end
        end

        
        %% ===================== 系统操作 =====================
        function clearSystemMatrices(obj)
            % 清空系统矩阵
            obj.GlobalStiffness = sparse(obj.NumTotalDofs, obj.NumTotalDofs);
            obj.GlobalMass = sparse(obj.NumTotalDofs, obj.NumTotalDofs);
        end
        
        %% ===================== 信息查询 =====================
        function info = get.SystemInfo(obj)
            % 获取系统基本信息
            info = struct();
            info.NumNodes = obj.NumTotalNodes;
            info.NumDofs = obj.NumTotalDofs;
            info.NumEntities = length(obj.Entities);
            info.NumMaterials = length(obj.Materials);
            
            if ~isempty(obj.GlobalStiffness)
                info.NonzeroStiffness = nnz(obj.GlobalStiffness);
                info.StiffnessSparsity = info.NonzeroStiffness / (info.NumDofs^2);
            end
        end
        
        function displaySystemInfo(obj)
            % 显示系统信息
            info = obj.SystemInfo;
            fprintf('=== 有限元系统信息 ===\n');
            fprintf('节点数: %d\n', info.NumNodes);
            fprintf('自由度数: %d\n', info.NumDofs);
            fprintf('实体数: %d\n', info.NumEntities);
            fprintf('材料数: %d\n', info.NumMaterials);
            
            if isfield(info, 'NonzeroStiffness')
                fprintf('刚度矩阵非零元: %d (稀疏度: %.4f%%)\n', ...
                    info.NonzeroStiffness, info.StiffnessSparsity * 100);
            end
        end
    end
end